
<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta charset="utf-8" />
    <title>smp &#8212; Texas A&amp;M Oil spill / Outfall Calculator 2.1.0 documentation</title>
    <link rel="stylesheet" href="../_static/alabaster.css" type="text/css" />
    <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
    <script id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
    <script src="../_static/jquery.js"></script>
    <script src="../_static/underscore.js"></script>
    <script src="../_static/doctools.js"></script>
    <script src="../_static/language_data.js"></script>
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
   
  <link rel="stylesheet" href="../_static/custom.css" type="text/css" />
  
  
  <meta name="viewport" content="width=device-width, initial-scale=0.9, maximum-scale=0.9" />

  </head><body>
  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          

          <div class="body" role="main">
            
  <h1>Source code for smp</h1><div class="highlight"><pre>
<span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Stratified Multiphase Plume</span>
<span class="sd">===========================</span>

<span class="sd">This module contains the numerical solution for the `stratified_plume_model`</span>
<span class="sd">module.  Some of the general tools for handling the multiphase components are</span>
<span class="sd">contained in `dispersed_phases`.</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="c1"># S. Socolofsky, November 2014, Texas A&amp;M University &lt;socolofs@tamu.edu&gt;.</span>

<span class="kn">from</span> <span class="nn">__future__</span> <span class="kn">import</span> <span class="p">(</span><span class="n">absolute_import</span><span class="p">,</span> <span class="n">division</span><span class="p">,</span> <span class="n">print_function</span><span class="p">)</span>

<span class="kn">from</span> <span class="nn">tamoc</span> <span class="kn">import</span> <span class="n">seawater</span>
<span class="kn">from</span> <span class="nn">tamoc</span> <span class="kn">import</span> <span class="n">dispersed_phases</span>

<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">scipy</span> <span class="kn">import</span> <span class="n">integrate</span>
<span class="kn">from</span> <span class="nn">scipy.optimize</span> <span class="kn">import</span> <span class="n">fsolve</span>

<span class="c1"># ----------------------------------------------------------------------------</span>
<span class="c1"># Derivatives of the system of ODEs</span>
<span class="c1"># ----------------------------------------------------------------------------</span>

<div class="viewcode-block" id="derivs_inner"><a class="viewcode-back" href="../autodoc/spm/smp.derivs_inner.html#smp.derivs_inner">[docs]</a><span class="k">def</span> <span class="nf">derivs_inner</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">yi</span><span class="p">,</span> <span class="n">yo</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculate the derivatives for the system of ODEs for the inner plume</span>
<span class="sd">    </span>
<span class="sd">    Calculates the right-hand-side of the system of ODEs for the inner plume</span>
<span class="sd">    state space.  These equations follow Socolofsky et al. (2008) very </span>
<span class="sd">    closely, with the exception that multiple dispersed phase particles are</span>
<span class="sd">    allowed within the inner plume.  Heat transfer between the dispersed</span>
<span class="sd">    and continuous phase is also added.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    z : float</span>
<span class="sd">        Current value for the independent variable (depth in m).</span>
<span class="sd">    y : ndarray</span>
<span class="sd">        Current value for the inner plume state space vector.</span>
<span class="sd">    yi : `InnerPlume`</span>
<span class="sd">        Object for manipulating the inner plume state space</span>
<span class="sd">    yo : `OuterPlume`</span>
<span class="sd">        Object for manipulating the outer plume state space</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase local</span>
<span class="sd">        conditions and behavior.</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    neighbor : `scipy.interpolate.interp1d` object</span>
<span class="sd">        Container holding the latest solution for the outer plume state</span>
<span class="sd">        space</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    yp : ndarray</span>
<span class="sd">        A vector of the derivatives of the inner plume state space.</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    stratified_plume_model.InnerPlume, stratified_plume_model.OuterPlume, </span>
<span class="sd">    stratified_plume_model.inner_main, calculate</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    It is important that the inner plume entrains fluid from either the </span>
<span class="sd">    ambient water (whenever the outer plume is not present) or the outer </span>
<span class="sd">    plume (whenever it is shrouding the inner plume).  This is accomplished</span>
<span class="sd">    in `stratified_plume_model.OuterPlume`:  if there is no outer plume </span>
<span class="sd">    segment, then the ambient conditions are stored in the outer plume </span>
<span class="sd">    variables.  Thus, `yo.c[i]` is equivalent to `ca[i]` when there is no </span>
<span class="sd">    outer plume.  This behavior is true for temperature, salinity, density</span>
<span class="sd">    and concentration.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Set up the output from the fuction to have the right size and type</span>
    <span class="n">yp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">yi</span><span class="o">.</span><span class="n">len</span><span class="p">,</span> <span class="mi">1</span><span class="p">))</span>
    
    <span class="c1"># Update the inner plume object with the corrent solution and compute</span>
    <span class="c1"># the inner plume shear entrainment coefficient</span>
    <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    
    <span class="c1"># Update the outer plume object at the current depth</span>
    <span class="k">if</span> <span class="n">z</span> <span class="o">&lt;</span> <span class="n">np</span><span class="o">.</span><span class="n">min</span><span class="p">(</span><span class="n">neighbor</span><span class="o">.</span><span class="n">x</span><span class="p">):</span>
        <span class="c1"># This plume is above any existing outer plumes</span>
        <span class="n">yo</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">len</span><span class="p">),</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Interpolate the outer plume solution to the current depth</span>
        <span class="n">yo</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">(</span><span class="n">z</span><span class="p">),</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="p">)</span>
    
    <span class="c1"># Conservation of mass</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> \
            <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span>
            
    <span class="c1"># Conservation of momentum</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">/</span> <span class="n">p</span><span class="o">.</span><span class="n">gamma_i</span> <span class="o">*</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">g</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> \
            <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">Fb</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">lambda_2</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">yi</span><span class="o">.</span><span class="n">Xi</span><span class="p">)</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">-</span> \
            <span class="n">yi</span><span class="o">.</span><span class="n">rho</span><span class="p">))</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> \
            <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">u</span><span class="p">)</span>
            
    <span class="c1"># Conservation of salinity</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> \
            <span class="n">yo</span><span class="o">.</span><span class="n">s</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">s</span><span class="p">)</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">s</span>
            
    <span class="c1"># Conservation of continuous phase fluid heat</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="p">(</span><span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> \
            <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">T</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> \
            <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
            
    <span class="c1"># Conservation equations for each dispersed phase</span>
    <span class="n">idx</span> <span class="o">=</span> <span class="mi">4</span>
    
    <span class="c1"># Track the mass dissolving into the continuous phase</span>
    <span class="n">delDiss</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">nchems</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">np</span><span class="p">):</span>
        <span class="n">delDiss_p</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">nchems</span><span class="p">)</span>
        
        <span class="k">if</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">issoluble</span><span class="p">:</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">nchems</span><span class="p">):</span>
                
                <span class="c1"># Conservation of particle mass for soluble particles </span>
                <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">A</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nb0</span> <span class="o">/</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> 
                            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">us</span><span class="p">)</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">beta</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">*</span> 
                            <span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">Cs</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">yi</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">]))</span>
                <span class="n">delDiss</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">+=</span> <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span>
                <span class="n">delDiss_p</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">+=</span> <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span>
                
                <span class="c1"># Update continuous phase temperature with heat of solution</span>
                <span class="n">yp</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">+=</span> <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">neg_dH_solR</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">Ru</span> <span class="o">/</span> \
                        <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">M</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
                <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            
        <span class="k">else</span><span class="p">:</span>
            <span class="c1"># Conservation of particle mass for insoluble particles</span>
            <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Conservation of particle heat including dissolution mass transfer</span>
        <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">A</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nb0</span> <span class="o">/</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> 
                  <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">us</span><span class="p">)</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">rho_p</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">cp</span> <span class="o">*</span> \
                  <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">beta_T</span> <span class="o">*</span> <span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span> <span class="o">-</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> \
                  <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">delDiss_p</span><span class="p">)</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">cp</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span>
        
        <span class="c1"># Take the heat leaving the particle and put it in the continuous </span>
        <span class="c1"># phase fluid</span>
        <span class="n">yp</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">-=</span> <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> 
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Track the age of each particle by following its advection</span>
        <span class="k">if</span> <span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">us</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
            <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">/</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">us</span><span class="p">)</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Track the location of each particle relative to the plume centerline</span>
        <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">:</span><span class="n">idx</span><span class="o">+</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">3</span>
    
    <span class="c1"># Conservation equations for the dissolved constituents.</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">nchems</span><span class="p">):</span>
        <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> \
                  <span class="n">yo</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> \
                  <span class="n">yi</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span> <span class="n">delDiss</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
    
    <span class="c1"># z is positive downward (depth)</span>
    <span class="k">return</span> <span class="o">-</span><span class="n">yp</span></div>

<div class="viewcode-block" id="derivs_outer"><a class="viewcode-back" href="../autodoc/spm/smp.derivs_outer.html#smp.derivs_outer">[docs]</a><span class="k">def</span> <span class="nf">derivs_outer</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">yi</span><span class="p">,</span> <span class="n">yo</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculate the derivatives for the system of ODEs for the outer plume</span>
<span class="sd">    </span>
<span class="sd">    Calculates the right-hand-side of the system of ODEs for the outer plume</span>
<span class="sd">    state space.  These equations follow those in Socolofsky et al. (2008).</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    z : float</span>
<span class="sd">        Current value for the independent variable (depth in m).</span>
<span class="sd">    y : ndarray</span>
<span class="sd">        Current value for the outer plume state space vector.</span>
<span class="sd">    yi : `InnerPlume`</span>
<span class="sd">        Object for manipulating the inner plume state space</span>
<span class="sd">    yo : `OuterPlume`</span>
<span class="sd">        Object for manipulating the outer plume state space</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase local</span>
<span class="sd">        conditions and behavior.</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    neighbor : `scipy.interpolate.interp1d` object</span>
<span class="sd">        Container holding the latest solution for the outer plume state</span>
<span class="sd">        space</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    yp : ndarray</span>
<span class="sd">        A vector of the derivatives of the outer plume state space.</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    stratified_plume_model.InnerPlume, stratified_plume_model.OuterPlume, </span>
<span class="sd">    stratified_plume_model.outer_main, calculate</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Set up the output from the function to have the correct size and type</span>
    <span class="n">yp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">yo</span><span class="o">.</span><span class="n">len</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
    
    <span class="c1"># Update the inner plume object at the current depth and the inner plume</span>
    <span class="c1"># shear entrainment coefficient</span>
    <span class="k">if</span> <span class="n">z</span> <span class="o">&gt;</span> <span class="n">np</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">neighbor</span><span class="o">.</span><span class="n">x</span><span class="p">):</span>
        <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">len</span><span class="p">),</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">(</span><span class="n">z</span><span class="p">),</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    
    <span class="c1"># Update the outer plume object with the current solution</span>
    <span class="n">yo</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="p">)</span>
    
    <span class="c1"># Conservation of Mass:</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> \
            <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_3</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span>
    
    <span class="c1"># Conservation of Momentum:</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.</span> <span class="o">/</span> <span class="n">p</span><span class="o">.</span><span class="n">gamma_o</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">g</span> <span class="o">*</span> <span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="o">-</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> \
            <span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">-</span> <span class="n">yo</span><span class="o">.</span><span class="n">rho</span><span class="p">)</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> \
            <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">u</span><span class="p">)</span>
    
    <span class="c1"># Conservation of Salinity:</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">s</span> <span class="o">+</span> \
            <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">s</span><span class="p">)</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_3</span> <span class="o">*</span> \
            <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">Sa</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">s</span>
    
    <span class="c1"># Conservation of Heat:</span>
    <span class="n">yp</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="p">(</span><span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> \
            <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">T</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span><span class="p">)</span> <span class="o">+</span> \
            <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_3</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">Ta</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
    
    <span class="c1"># Conservation of tracked chemical constituents:</span>
    <span class="n">idx</span> <span class="o">=</span> <span class="mi">4</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">nchems</span><span class="p">):</span>
        <span class="n">yp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span><span class="p">)</span> <span class="o">*</span> \
                <span class="n">yo</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> \
                <span class="n">p</span><span class="o">.</span><span class="n">alpha_3</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">u</span> <span class="o">*</span> <span class="n">yo</span><span class="o">.</span><span class="n">ca</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
    
    <span class="c1"># z is positive downward (depth)</span>
    <span class="k">return</span> <span class="n">yp</span></div>


<span class="c1"># ----------------------------------------------------------------------------</span>
<span class="c1"># Main integration controller</span>
<span class="c1"># ----------------------------------------------------------------------------</span>

<div class="viewcode-block" id="calculate"><a class="viewcode-back" href="../autodoc/spm/smp.calculate.html#smp.calculate">[docs]</a><span class="k">def</span> <span class="nf">calculate</span><span class="p">(</span><span class="n">yi</span><span class="p">,</span> <span class="n">yo</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">,</span> <span class="n">derivs</span><span class="p">,</span> <span class="n">z0</span><span class="p">,</span> <span class="n">y0</span><span class="p">,</span> <span class="n">zf</span><span class="p">,</span> 
              <span class="n">z_dir</span><span class="p">,</span> <span class="n">delta_z</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Integrate an inner or outer plume segment from `z0` to `zf`</span>
<span class="sd">    </span>
<span class="sd">    Integrate the inner or outer plume over the range from `z0` to `zf`, </span>
<span class="sd">    integrating in the direction (positive or negative) given by `z_dir`.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    yi : `InnerPlume`</span>
<span class="sd">        Object for manipulating the inner plume state space</span>
<span class="sd">    yo : `OuterPlume`</span>
<span class="sd">        Object for manipulating the outer plume state space</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase local</span>
<span class="sd">        conditions and behavior.</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    neighbor : `scipy.interpolate.interp1d` object</span>
<span class="sd">        Container holding the latest solution for the outer plume state</span>
<span class="sd">        space</span>
<span class="sd">    derivs : function handle</span>
<span class="sd">        Pointer to the function where the derivatives of the ODE system are</span>
<span class="sd">        stored.  Should be either `smp.derivs_inner` or `smp.derivs_outer`.</span>
<span class="sd">    z0 : float</span>
<span class="sd">        Initial depth (m)</span>
<span class="sd">    y0 : ndarray</span>
<span class="sd">        Initial values of the state space vector</span>
<span class="sd">    zf : float</span>
<span class="sd">        Final depth to calculate (m)</span>
<span class="sd">    z_dir : float</span>
<span class="sd">        Direction (+1 or -1) to integrate the vertical coordinate.  The inner</span>
<span class="sd">        plume integrates in the negative z-direction (to shallower depths),</span>
<span class="sd">        and the outer plume integrates in the positive z-direction (to </span>
<span class="sd">        greater depths).</span>
<span class="sd">    delta_z : float</span>
<span class="sd">        Maximum step size to use in the simulation (m).  The ODE solver </span>
<span class="sd">        in `calculate` is set up with adaptive step size integration, so </span>
<span class="sd">        this value determines the largest step size in the output data, but </span>
<span class="sd">        not the numerical stability of the calculation.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    z : ndarray</span>
<span class="sd">        Vector of elevations where the inner plume solution is obtained (m).</span>
<span class="sd">    y : ndarray</span>
<span class="sd">        Matrix of inner plume state space solutions.  Each row corresponds to</span>
<span class="sd">        a depth in `z`.</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    derivs_inner, derivs_outer, stratified_plume_model.Model, </span>
<span class="sd">    stratified_plume_model.inner_main, stratified_plume_model.outer_main, </span>
<span class="sd">    stratified_plume_model.InnerPlume, stratified_plume_model.OuterPlume</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Create the integrator object:  use &quot;vode&quot; with &quot;backward </span>
    <span class="c1"># differentiation formula&quot; for stiff ODEs</span>
    <span class="n">r</span> <span class="o">=</span> <span class="n">integrate</span><span class="o">.</span><span class="n">ode</span><span class="p">(</span><span class="n">derivs</span><span class="p">)</span><span class="o">.</span><span class="n">set_integrator</span><span class="p">(</span><span class="s1">&#39;vode&#39;</span><span class="p">,</span> <span class="n">method</span><span class="o">=</span><span class="s1">&#39;bdf&#39;</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="mf">1.e-6</span><span class="p">,</span> 
        <span class="n">rtol</span><span class="o">=</span><span class="mf">1e-3</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="mi">5</span><span class="p">,</span> <span class="n">max_step</span><span class="o">=</span><span class="n">delta_z</span><span class="p">)</span>
    
    <span class="c1"># Initialize the state space</span>
    <span class="n">r</span><span class="o">.</span><span class="n">set_initial_value</span><span class="p">(</span><span class="n">y0</span><span class="p">,</span> <span class="n">z0</span><span class="p">)</span>
    
    <span class="c1"># Set passing variables for derivs method</span>
    <span class="n">r</span><span class="o">.</span><span class="n">set_f_params</span><span class="p">(</span><span class="n">yi</span><span class="p">,</span> <span class="n">yo</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">)</span>
    
    <span class="c1"># Create vectors (using the list data type) to store the solution</span>
    <span class="n">z</span> <span class="o">=</span> <span class="p">[</span><span class="n">z0</span><span class="p">]</span>
    <span class="n">y</span> <span class="o">=</span> <span class="p">[</span><span class="n">y0</span><span class="p">]</span>
    
    <span class="c1"># Integrate to zf unless the solution stops naturally earlier</span>
    <span class="n">k</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">psteps</span> <span class="o">=</span> <span class="mi">30</span>
    <span class="n">stop</span> <span class="o">=</span> <span class="kc">False</span>
    <span class="k">while</span> <span class="n">r</span><span class="o">.</span><span class="n">successful</span><span class="p">()</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">stop</span><span class="p">:</span>
        
        <span class="c1"># Print progress to the screen</span>
        <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">remainder</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">float</span><span class="p">(</span><span class="n">k</span><span class="p">),</span> <span class="n">psteps</span><span class="p">)</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
            <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;    Depth:  </span><span class="si">%g</span><span class="s1"> (m), k: </span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">z</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">k</span><span class="p">))</span>
        
        <span class="c1"># Perform one step of the integration</span>
        <span class="n">r</span><span class="o">.</span><span class="n">integrate</span><span class="p">(</span><span class="n">zf</span><span class="p">,</span> <span class="n">step</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
        
        <span class="c1"># Store the results</span>
        <span class="k">if</span> <span class="n">derivs</span> <span class="o">==</span> <span class="n">derivs_inner</span><span class="p">:</span>
            <span class="c1"># Store the correct temperature for the particles after heat </span>
            <span class="c1"># transfer turns off</span>
            <span class="n">r</span> <span class="o">=</span> <span class="n">correct_temperature</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">yi</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
        <span class="n">z</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">t</span><span class="p">)</span>
        <span class="n">y</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">)</span>
        <span class="n">k</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Evaluate the stop criteria</span>
        <span class="k">if</span> <span class="n">r</span><span class="o">.</span><span class="n">successful</span><span class="p">():</span>
            <span class="c1"># Check if we reached the free surface</span>
            <span class="k">if</span> <span class="n">z</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">*</span> <span class="n">z_dir</span> <span class="o">&gt;=</span> <span class="n">zf</span> <span class="o">*</span> <span class="n">z_dir</span><span class="p">:</span>
                <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
            <span class="c1"># Check if the momentum went negative</span>
            <span class="k">if</span> <span class="n">y</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span> <span class="o">&lt;</span> <span class="mf">0.</span><span class="p">:</span>
                <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
            <span class="c1"># Check if the progress stopped</span>
            <span class="k">if</span> <span class="n">z</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">==</span> <span class="n">z</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]:</span>
                <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
    
    <span class="c1"># Convert solution to numpy arrays.</span>
    <span class="n">z</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">z</span><span class="p">)</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">y</span><span class="p">)</span>
    
    <span class="c1"># Remove any part of the solution with a negative momentum</span>
    <span class="n">rows</span> <span class="o">=</span> <span class="n">y</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;=</span> <span class="mi">0</span>
    <span class="n">z</span> <span class="o">=</span> <span class="n">z</span><span class="p">[</span><span class="n">rows</span><span class="p">]</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">y</span><span class="p">[</span><span class="n">rows</span><span class="p">,:]</span>
    
    <span class="c1"># Return the solution</span>
    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">remainder</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">float</span><span class="p">(</span><span class="n">k</span><span class="p">),</span> <span class="n">psteps</span><span class="p">)</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
        <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;    Depth:  </span><span class="si">%g</span><span class="s1"> (m), k: </span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">z</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">k</span><span class="p">))</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">y</span><span class="p">)</span></div>

<div class="viewcode-block" id="correct_temperature"><a class="viewcode-back" href="../autodoc/spm/smp.correct_temperature.html#smp.correct_temperature">[docs]</a><span class="k">def</span> <span class="nf">correct_temperature</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">yi</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Make sure the correct temperature is stored in the state space solution</span>
<span class="sd">    </span>
<span class="sd">    When the dispersed phase particles equilibrate to their surrounding </span>
<span class="sd">    temperature, heat transfer is turned off by the methods in </span>
<span class="sd">    `dispersed_phases.Particle`.  This is needed to prevent numerical</span>
<span class="sd">    oscillations as the particles become small.  Unfortunately, it is not as</span>
<span class="sd">    easy to make the numerical solution output the correct result once </span>
<span class="sd">    particle temperature effectively stops being a state space variable.  </span>
<span class="sd">    </span>
<span class="sd">    Once heat transfer is turned off, all of the model methods use the </span>
<span class="sd">    correct temperature (e.g., the ambient temperature) in all of the </span>
<span class="sd">    equations coupled to the heat transfer equation and in all equations </span>
<span class="sd">    involving particle temperature.  </span>
<span class="sd">    </span>
<span class="sd">    In order to prevent the state space variable for particle temperature </span>
<span class="sd">    from blowing up as the mass goes to zero, we also continue to adjust the</span>
<span class="sd">    particle heat in the ODE solution to maintain a constant temperature. </span>
<span class="sd">    This is done by setting `beta_T = 0`.  This is merely a numerical trick, </span>
<span class="sd">    as all equations using the particle temperature know to use the ambient</span>
<span class="sd">    temperature when this is the case.  </span>
<span class="sd">    </span>
<span class="sd">    Hence, the purpose of this function is to simply overwrite the state </span>
<span class="sd">    space solution containing the particle heat (returned by the ODE solver</span>
<span class="sd">    to maintain a constant particle temperature) with the correct particle</span>
<span class="sd">    heat yielding the ambient temperature for the particle temperature.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    r : `scipy.integrate.ode` object</span>
<span class="sd">        ODE solution containing the currect values of the state space (e.g., </span>
<span class="sd">        `r.y`).</span>
<span class="sd">    yi : `InnerPlume`</span>
<span class="sd">        Object for manipulating the inner plume state space</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase local</span>
<span class="sd">        conditions and behavior.</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    r : `scipy.integrate.ode` object</span>
<span class="sd">        Returns the original ODE object with the corrected solution stored</span>
<span class="sd">        in the public x and y.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Update the inner plume state space with the current solution.  This </span>
    <span class="c1"># will set the correct particle temperature in the attributes </span>
    <span class="c1"># yi.particles[].T.  If heat transfer is still turned on, the answer will</span>
    <span class="c1"># be the value computed from the state space r.x and r.y; if heat </span>
    <span class="c1"># transfer is turned off, the answer will be the ambient fluid</span>
    <span class="c1"># temperature (e.g., Ti).</span>
    <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">t</span><span class="p">,</span> <span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    
    <span class="c1"># Find the heat conservation equation in the inner plume state space and</span>
    <span class="c1"># replace the heat with the correct value so that r.y always yields the</span>
    <span class="c1"># particle temperature determined above.</span>
    <span class="n">idx</span> <span class="o">=</span> <span class="mi">4</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">nc</span>
        <span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">m</span><span class="p">)</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nb0</span> <span class="o">*</span> \
                       <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">cp</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span>
        <span class="c1"># Advance for heat, time, and position</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span> <span class="o">+</span> <span class="mi">1</span> <span class="o">+</span> <span class="mi">3</span>
    
    <span class="c1"># Return the corrected solution</span>
    <span class="k">return</span> <span class="n">r</span></div>

<span class="c1"># ----------------------------------------------------------------------------</span>
<span class="c1"># General tools used throughout the model simulation</span>
<span class="c1"># ----------------------------------------------------------------------------</span>

<div class="viewcode-block" id="cp_model"><a class="viewcode-back" href="../autodoc/spm/smp.cp_model.html#smp.cp_model">[docs]</a><span class="k">def</span> <span class="nf">cp_model</span><span class="p">(</span><span class="n">epsilon</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">rho_a</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">rho_r</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">u</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Continuous peeling model of Crounse (2000)</span>
<span class="sd">    </span>
<span class="sd">    Computes the local peeling flux from the continuous model of Crounse</span>
<span class="sd">    (2000)</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    epsilon : float</span>
<span class="sd">        Continuous peeling model calibration factor (--).</span>
<span class="sd">    particles : list of `single_particle_model.Particle` objects</span>
<span class="sd">        Iterable list of dbm class objects describing each dispersed phase.</span>
<span class="sd">    rho_a : float</span>
<span class="sd">        Local density of ambient fluid outside plume (kg/m^3).</span>
<span class="sd">    rho : float</span>
<span class="sd">        Local density of plume fluid (kg/m^3)</span>
<span class="sd">    g : float</span>
<span class="sd">        Acceleration of gravity (m/s^2).</span>
<span class="sd">    rho_r : float</span>
<span class="sd">        Model reference density (kg/m^3).</span>
<span class="sd">    b : float</span>
<span class="sd">        Local plume half-width (m)</span>
<span class="sd">    u : float</span>
<span class="sd">        Local plume fluid velocity (m/s)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute a buoyancy flux weighted average of the slip velocity</span>
    <span class="n">us</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">us</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">))])</span>
    <span class="n">us</span> <span class="o">=</span> <span class="n">dispersed_phases</span><span class="o">.</span><span class="n">bf_average</span><span class="p">(</span><span class="n">particles</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">rho_r</span><span class="p">,</span> <span class="n">us</span><span class="p">)</span>
    
    <span class="c1"># Return the peeling flux</span>
    <span class="k">return</span> <span class="n">epsilon</span> <span class="o">*</span> <span class="p">(</span><span class="n">us</span> <span class="o">/</span> <span class="n">u</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">g</span> <span class="o">*</span> <span class="p">(</span><span class="n">rho_a</span> <span class="o">-</span> <span class="n">rho</span><span class="p">)</span> <span class="o">/</span> <span class="n">rho_r</span> <span class="o">*</span> \
           <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">u</span></div>

<span class="c1"># ----------------------------------------------------------------------------</span>
<span class="c1"># Functions to compute inner plume initial conditions</span>
<span class="c1"># ---------------------------------------------------------------------------</span>

<div class="viewcode-block" id="main_ic"><a class="viewcode-back" href="../autodoc/spm/smp.main_ic.html#smp.main_ic">[docs]</a><span class="k">def</span> <span class="nf">main_ic</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">z0</span><span class="p">,</span> <span class="n">R</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the initial conditions for the inner plume state space</span>
<span class="sd">    </span>
<span class="sd">    Compute the initial conditions at the release location for the inner</span>
<span class="sd">    plume state space</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the single bubble model simulation.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase initial</span>
<span class="sd">        conditions</span>
<span class="sd">    p : `stratified_plume_model.ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    z0 : float</span>
<span class="sd">        Depth of the release point (m)</span>
<span class="sd">    R : float</span>
<span class="sd">        Radius of the release port (m)</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    z : float</span>
<span class="sd">        Initial depth for the simulation (m)</span>
<span class="sd">    q : ndarray</span>
<span class="sd">        Initial value of the inner plume state space</span>
<span class="sd">    chem_names : str list</span>
<span class="sd">        List of the chemicals in the dispersed phase composition that are </span>
<span class="sd">        undergoing dissolution</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the initial volume flux at the source</span>
    <span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">z</span><span class="p">,</span> <span class="n">Ta</span><span class="p">,</span> <span class="n">Sa</span><span class="p">,</span> <span class="n">P</span><span class="p">,</span> <span class="n">rho</span> <span class="o">=</span> <span class="n">dispersed_phases</span><span class="o">.</span><span class="n">zfe_volume_flux</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> 
                               <span class="n">particles</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">z0</span><span class="p">,</span> <span class="n">R</span><span class="p">)</span>
    
    <span class="c1"># Get the dispersed phase chemical components</span>
    <span class="n">chem_names</span> <span class="o">=</span> <span class="n">dispersed_phases</span><span class="o">.</span><span class="n">get_chem_names</span><span class="p">(</span><span class="n">particles</span><span class="p">)</span>
        
    <span class="c1"># Build the initial state space with these initial values</span>
    <span class="n">z</span><span class="p">,</span> <span class="n">y</span> <span class="o">=</span> <span class="n">inner_plume_ic</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">z</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">Sa</span><span class="p">,</span> <span class="n">Ta</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">)</span>
    
    <span class="c1"># Return the initial depth, state space, and list of chem_names</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">)</span></div>


<div class="viewcode-block" id="inner_plume_ic"><a class="viewcode-back" href="../autodoc/spm/smp.inner_plume_ic.html#smp.inner_plume_ic">[docs]</a><span class="k">def</span> <span class="nf">inner_plume_ic</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">z</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">S</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Build the inner plume state space given the initial conditions</span>
<span class="sd">    </span>
<span class="sd">    Constructs the state space for the inner plume from the initial values for</span>
<span class="sd">    Q, J, concentrations, and particle properties.  The state space vector is</span>
<span class="sd">    organized as follows:</span>
<span class="sd">    </span>
<span class="sd">        y[0] = Q : Flow rate of entrained fluid</span>
<span class="sd">        y[1] = J : Momentum flux of entrained fluid</span>
<span class="sd">        y[2] = S : Salinity flux of entrained fluid</span>
<span class="sd">        y[3] = H : Heat flux of entrained fluid</span>
<span class="sd">        y[4:4 + np * (nchems + 1)] : Dispersed phase mass and heat fluxes</span>
<span class="sd">        y[5 + np * (nchems + 1):] : Mass fluxes of the dissolved components</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the single bubble model simulation.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase initial</span>
<span class="sd">        conditions</span>
<span class="sd">    p : `stratified_plume_model.ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    z : float</span>
<span class="sd">        Depth of the release point (m)</span>
<span class="sd">    Q : float</span>
<span class="sd">        Initial volume flux of entrained seawater (m^3/s)</span>
<span class="sd">    A : float</span>
<span class="sd">        Cross-sectional area of the discharge (m^2)</span>
<span class="sd">    S : float</span>
<span class="sd">        Salinity of the entrained seawater (psu)</span>
<span class="sd">    T : float</span>
<span class="sd">        Temperature of the entrained seawater (K)</span>
<span class="sd">    chem_names : string list</span>
<span class="sd">        List of the names of the chemicals that will be tracked in the </span>
<span class="sd">        dissolved phase</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    z : float</span>
<span class="sd">        Depth at the initial point of the plume (m)</span>
<span class="sd">    y : ndarray</span>
<span class="sd">        Initial inner plume state space (see description above)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Sequentially build the inner plume state space</span>
    <span class="n">y</span> <span class="o">=</span> <span class="p">[</span><span class="n">Q</span><span class="p">,</span> <span class="n">Q</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">A</span><span class="p">,</span> <span class="n">S</span> <span class="o">*</span> <span class="n">Q</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="n">T</span> <span class="o">*</span> <span class="n">Q</span><span class="p">]</span>
    
    <span class="c1"># Add in the state space of the multiphase components</span>
    <span class="n">nb0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">))</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>
        <span class="n">nb0</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nb0</span>
    <span class="n">y</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">dispersed_phases</span><span class="o">.</span><span class="n">particles_state_space</span><span class="p">(</span><span class="n">particles</span><span class="p">,</span> <span class="n">nb0</span><span class="p">))</span>
    
    <span class="c1"># And the mass fluxes of dissolved components</span>
    <span class="n">ca</span> <span class="o">=</span> <span class="n">profile</span><span class="o">.</span><span class="n">get_values</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">)</span>
    <span class="n">y</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">ca</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    
    <span class="c1"># Return the initial state space</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">y</span><span class="p">))</span></div>


<span class="c1"># ----------------------------------------------------------------------------</span>
<span class="c1"># Functions to compute outer plume initial conditions</span>
<span class="c1"># ---------------------------------------------------------------------------</span>

<div class="viewcode-block" id="outer_surf"><a class="viewcode-back" href="../autodoc/spm/smp.outer_surf.html#smp.outer_surf">[docs]</a><span class="k">def</span> <span class="nf">outer_surf</span><span class="p">(</span><span class="n">yi</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the initial condition for the outer plume at the sea surface</span>
<span class="sd">    </span>
<span class="sd">    Computes the initial conditions for the first outer plume segment after</span>
<span class="sd">    the inner plume impinges on the free surface of the water body.  It is </span>
<span class="sd">    assumed that the inner plume had significant volume flux and that this </span>
<span class="sd">    first outer plume segment will be viable.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    yi : `stratified_plume_model.InnerPlume` object</span>
<span class="sd">        Object for manipulating the inner plume state space</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    z0 : float</span>
<span class="sd">        Initial depth of the outer plume segment (m).</span>
<span class="sd">    y0 : ndarray</span>
<span class="sd">        Initial dependent variables state space for the outer plume segment.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># The outer plume is a mixture of inner plume fluid and ambient fluid</span>
    <span class="c1"># entrained from the water surface</span>
    <span class="n">Q</span> <span class="o">=</span> <span class="p">(</span><span class="mf">1.</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">fe</span><span class="p">)</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">Q</span>
    <span class="n">T</span> <span class="o">=</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">T</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ta</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">fe</span><span class="p">)</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">Q</span> <span class="o">/</span> <span class="n">Q</span>
    <span class="n">s</span> <span class="o">=</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">s</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Sa</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">fe</span><span class="p">)</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">Q</span> <span class="o">/</span> <span class="n">Q</span>
    <span class="n">c</span> <span class="o">=</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">c</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">ca</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">fe</span><span class="p">)</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">Q</span> <span class="o">/</span> <span class="n">Q</span>
    <span class="n">rho</span> <span class="o">=</span> <span class="n">seawater</span><span class="o">.</span><span class="n">density</span><span class="p">(</span><span class="n">T</span><span class="p">,</span> <span class="n">s</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">P</span><span class="p">)</span>
    
    <span class="c1"># Use a Froude number approach to set the initial width and velocity</span>
    <span class="n">u</span> <span class="o">=</span> <span class="n">outer_fr</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">rho_a</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">g</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">Fro_0</span><span class="p">)</span>
    
    <span class="c1"># Calculate the outer plume state space variables</span>
    <span class="n">y0</span> <span class="o">=</span> <span class="p">[]</span>
    <span class="n">Q</span> <span class="o">=</span> <span class="o">-</span><span class="n">Q</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Q</span><span class="p">)</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Q</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">))</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">s</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="n">T</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">c</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    
    <span class="c1"># Return the outer plume initial condition</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">z</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">y0</span><span class="p">))</span></div>

<div class="viewcode-block" id="outer_dis"><a class="viewcode-back" href="../autodoc/spm/smp.outer_dis.html#smp.outer_dis">[docs]</a><span class="k">def</span> <span class="nf">outer_dis</span><span class="p">(</span><span class="n">yi</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">,</span> <span class="n">z_0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the initial condition for the outer plume at the DMPR</span>
<span class="sd">    </span>
<span class="sd">    Computes the initial conditions for the an outer plume segment at the </span>
<span class="sd">    depth of maximum plume rise (DMPR) following full dissolution of the </span>
<span class="sd">    dispersed phases.  </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    yi : `stratified_plume_model.InnerPlume` object</span>
<span class="sd">        Object for manipulating the inner plume state space.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase local</span>
<span class="sd">        conditions and behavior.</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    neighbor : `scipy.interpolate.interp1d` object</span>
<span class="sd">        Container holding the latest solution for the inner plume state</span>
<span class="sd">        space.</span>
<span class="sd">    z_0 : float</span>
<span class="sd">        Top of the inner plume calculation (m).</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    z0 : float</span>
<span class="sd">        Initial depth of the outer plume segment (m).</span>
<span class="sd">    y0 : ndarray</span>
<span class="sd">        Initial dependent variables state space for the outer plume segment.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Search for the maximum flux near the top of the plume</span>
    <span class="n">Qmax</span> <span class="o">=</span> <span class="n">neighbor</span><span class="o">.</span><span class="n">y</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">imax</span> <span class="o">=</span> <span class="mi">1</span>
    <span class="k">while</span> <span class="n">Qmax</span> <span class="o">&lt;</span> <span class="n">neighbor</span><span class="o">.</span><span class="n">y</span><span class="p">[</span><span class="n">imax</span><span class="p">,</span><span class="mi">0</span><span class="p">]:</span>
        <span class="n">Qmax</span> <span class="o">=</span> <span class="n">neighbor</span><span class="o">.</span><span class="n">y</span><span class="p">[</span><span class="n">imax</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">imax</span> <span class="o">+=</span> <span class="mi">1</span>
    
    <span class="c1"># Since most of this fluid will be regained as the outer plume descends</span>
    <span class="c1"># through Ep, take the initial volume flux as a small fraction (given</span>
    <span class="c1"># by model parameter qdis_ic)</span>
    <span class="n">Q</span> <span class="o">=</span> <span class="n">p</span><span class="o">.</span><span class="n">qdis_ic</span> <span class="o">*</span> <span class="n">Qmax</span>
    
    <span class="c1"># Get the local plume properties at the top of the plume</span>
    <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z_0</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">(</span><span class="n">z_0</span><span class="p">),</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    <span class="n">rho</span> <span class="o">=</span> <span class="n">yi</span><span class="o">.</span><span class="n">rho</span>
    
    <span class="c1"># Use a Froude number approach to set the initial width and velocity</span>
    <span class="n">u</span> <span class="o">=</span> <span class="n">outer_fr</span><span class="p">(</span><span class="mf">0.05</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">rho_a</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">g</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">Fro_0</span><span class="p">)</span>
    
    <span class="c1"># Calculate the outer plume state space variables</span>
    <span class="n">y0</span> <span class="o">=</span> <span class="p">[]</span>
    <span class="n">Q</span> <span class="o">=</span> <span class="o">-</span><span class="n">Q</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Q</span><span class="p">)</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Q</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">))</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">s</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    <span class="n">y0</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">c</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    
    <span class="c1"># Return the outer plume initial condition</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">z</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">y0</span><span class="p">))</span></div>

<div class="viewcode-block" id="outer_cpic"><a class="viewcode-back" href="../autodoc/spm/smp.outer_cpic.html#smp.outer_cpic">[docs]</a><span class="k">def</span> <span class="nf">outer_cpic</span><span class="p">(</span><span class="n">yi</span><span class="p">,</span> <span class="n">yo</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">,</span> <span class="n">z_0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the initial condition for the outer plume at depth</span>
<span class="sd">    </span>
<span class="sd">    Computes the initial conditions for the an outer plume segment within the </span>
<span class="sd">    reservoir body.  Part of the calculation determines whether or not the </span>
<span class="sd">    computed initial condition has enough downward momentum to be viable as </span>
<span class="sd">    an initial condition (e.g., whether or not it will be overwhelmed by the</span>
<span class="sd">    upward drag of the inner plume).</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    yi : `stratified_plume_model.InnerPlume` object</span>
<span class="sd">        Object for manipulating the inner plume state space.</span>
<span class="sd">    yo : `stratified_plume_model.OuterPlume` object</span>
<span class="sd">        Object for manipulating the outer plume state space.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `Particle` objects containing the dispersed phase local</span>
<span class="sd">        conditions and behavior.</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    neighbor : `scipy.interpolate.interp1d` object</span>
<span class="sd">        Container holding the latest solution for the inner plume state</span>
<span class="sd">        space.</span>
<span class="sd">    z_0 : float</span>
<span class="sd">        Top of the inner plume calculation (m).</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    z0 : float</span>
<span class="sd">        Initial depth of the outer plume segment (m).</span>
<span class="sd">    y0 : ndarray</span>
<span class="sd">        Initial dependent variables state space for the outer plume segment.</span>
<span class="sd">    flag : bool</span>
<span class="sd">        Outer plume viability flag:  `True` means the outer plume segment is</span>
<span class="sd">        viable and should be integrated; `False` means the outer plume </span>
<span class="sd">        segment is too weak and should be discarded, moving down the inner </span>
<span class="sd">        plume to calculate the next outer plume initial condition.</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    The iteration required to find a viable outer plume segment is conducted </span>
<span class="sd">    by the `stratified_plume_model.outer_main` function.  This function </span>
<span class="sd">    computes the initial conditions for one attempt to find an outer plume</span>
<span class="sd">    segment and reports back (through `flag`) on the success.</span>
<span class="sd">    </span>
<span class="sd">    There is one caveat to the above statement.  The model parameter </span>
<span class="sd">    `p.nwidths` determines the vertical scale over which this function may</span>
<span class="sd">    integrate to find the start to an outer plume, given as a integer number</span>
<span class="sd">    of times of the inner plume half-width.  This function starts by searching</span>
<span class="sd">    one half-width.  If `p.nwidths` is greater than one, it will continue to</span>
<span class="sd">    expand the search region.  The physical interpretation of `p.nwidths` is</span>
<span class="sd">    to set a reasonable upper bound on the diameter of eddies shed from the</span>
<span class="sd">    inner plume in the peeling region into the outer plume.  While the </span>
<span class="sd">    integral model does not have &quot;eddies&quot; per se, the search window size </span>
<span class="sd">    should still be representative of this type of length scale.  </span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Start the iteration counters</span>
    <span class="nb">iter</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">done</span> <span class="o">=</span> <span class="kc">False</span>
    
    <span class="c1"># Compute the outer plume initial conditions until the outer plume is </span>
    <span class="c1"># viable or until the maximum number of widths is integrated</span>
    <span class="k">while</span> <span class="ow">not</span> <span class="n">done</span> <span class="ow">and</span> <span class="nb">iter</span> <span class="o">&lt;</span> <span class="n">p</span><span class="o">.</span><span class="n">nwidths</span><span class="p">:</span>
        
        <span class="c1"># Update iteration counter</span>
        <span class="nb">iter</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Get the inner plume properties at the top of this peeling region</span>
        <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">z_0</span><span class="p">,</span> <span class="n">neighbor</span><span class="p">(</span><span class="n">z_0</span><span class="p">),</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
        
        <span class="c1"># Set the range to integrate to get the current peeling flux</span>
        <span class="n">z_upper</span> <span class="o">=</span> <span class="n">z_0</span>
        <span class="n">z_lower</span> <span class="o">=</span> <span class="n">z_0</span> <span class="o">+</span> <span class="nb">iter</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span>
        
        <span class="c1"># Check if the bottom of the reservoir is encountered.</span>
        <span class="k">if</span> <span class="n">z_lower</span> <span class="o">&gt;</span> <span class="n">profile</span><span class="o">.</span><span class="n">z_max</span><span class="p">:</span>
            <span class="n">z_lower</span> <span class="o">=</span> <span class="n">profile</span><span class="o">.</span><span class="n">z_max</span>
        
        <span class="c1"># Find the indices in the raw data for the inner plume solution close</span>
        <span class="c1"># to where z_upper and z_lower occur</span>
        <span class="n">i_upper</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">min</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="n">neighbor</span><span class="o">.</span><span class="n">x</span> <span class="o">&gt;=</span> <span class="n">z_upper</span><span class="p">)[</span><span class="mi">0</span><span class="p">])</span>
        <span class="n">i_lower</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="n">neighbor</span><span class="o">.</span><span class="n">x</span> <span class="o">&lt;=</span> <span class="n">z_lower</span><span class="p">)[</span><span class="mi">0</span><span class="p">])</span>
        
        <span class="c1"># Get the grid of elevations where we will integrate the solution to </span>
        <span class="c1"># obtain the initial flux for the outer plume.  This is needed </span>
        <span class="c1"># because the solution is so stiff:  if we integrated over a fixed</span>
        <span class="c1"># step size, we could easily miss dramatic changes in the solution.</span>
        <span class="c1"># Hence, we integrate over the steps in the numerical solution </span>
        <span class="c1"># itself.</span>
        <span class="n">n_grid</span> <span class="o">=</span> <span class="n">i_lower</span> <span class="o">-</span> <span class="n">i_upper</span> <span class="o">+</span> <span class="mi">3</span>
        <span class="n">zi</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n_grid</span><span class="p">)</span>
        <span class="n">zi</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">z_upper</span>
        <span class="n">zi</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">z_lower</span>
        <span class="n">zi</span><span class="p">[</span><span class="mi">1</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">neighbor</span><span class="o">.</span><span class="n">x</span><span class="p">[</span><span class="n">i_upper</span><span class="p">:</span><span class="n">i_lower</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span>
        
        <span class="c1"># Integrate the peeling fluid over this grid to get the total </span>
        <span class="c1"># contributions going into the outer plume</span>
        <span class="n">Q</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="n">tracer_vars</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="mi">2</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">nchems</span><span class="p">)</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">zi</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
            <span class="n">yi</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">zi</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="n">neighbor</span><span class="p">(</span><span class="n">zi</span><span class="p">[</span><span class="n">i</span><span class="p">]),</span> <span class="n">particles</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
            <span class="n">dz</span> <span class="o">=</span> <span class="n">zi</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">zi</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            <span class="n">Q</span> <span class="o">=</span> <span class="n">Q</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">dz</span>
            <span class="n">tracer_vars</span> <span class="o">=</span> <span class="n">tracer_vars</span> <span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">hstack</span><span class="p">((</span><span class="n">yi</span><span class="o">.</span><span class="n">s</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">T</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> \
                          <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">(),</span> <span class="n">yi</span><span class="o">.</span><span class="n">c</span><span class="p">))</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span> <span class="o">*</span> <span class="n">dz</span>
        
        <span class="c1"># Get the initial velocity of the peeling fluid using the modified </span>
        <span class="c1"># outer plume Froude number condition</span>
        <span class="n">T</span> <span class="o">=</span> <span class="n">tracer_vars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">/</span> <span class="p">(</span><span class="n">Q</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">())</span>
        <span class="n">s</span> <span class="o">=</span> <span class="n">tracer_vars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">/</span> <span class="n">Q</span>
        <span class="n">c</span> <span class="o">=</span> <span class="n">tracer_vars</span><span class="p">[</span><span class="mi">2</span><span class="p">:]</span> <span class="o">/</span> <span class="n">Q</span>
        <span class="n">rho</span> <span class="o">=</span> <span class="n">seawater</span><span class="o">.</span><span class="n">density</span><span class="p">(</span><span class="n">T</span><span class="p">,</span> <span class="n">s</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">P</span><span class="p">)</span>
        <span class="n">u</span> <span class="o">=</span> <span class="n">outer_fr</span><span class="p">(</span><span class="mf">0.05</span><span class="p">,</span> <span class="o">-</span><span class="n">Q</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="n">yi</span><span class="o">.</span><span class="n">rho_a</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">g</span><span class="p">,</span> <span class="n">p</span><span class="o">.</span><span class="n">Fro_0</span><span class="p">)</span>
        <span class="n">b</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">Q</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">Q</span><span class="p">)</span> <span class="o">*</span> <span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        <span class="n">dQdz</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">yi</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">p</span><span class="o">.</span><span class="n">alpha_1</span> <span class="o">*</span> <span class="p">(</span><span class="n">yi</span><span class="o">.</span><span class="n">u</span> <span class="o">+</span> <span class="n">p</span><span class="o">.</span><span class="n">c1</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">))</span> <span class="o">+</span> \
                <span class="n">p</span><span class="o">.</span><span class="n">alpha_2</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">))</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">b</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">alpha_3</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">yi</span><span class="o">.</span><span class="n">Ep</span>
        
        <span class="c1"># Check whether this outer plume segment will be viable</span>
        <span class="k">if</span> <span class="n">dQdz</span> <span class="o">&gt;</span> <span class="mi">0</span> <span class="ow">or</span> <span class="n">Q</span> <span class="o">&gt;</span> <span class="mi">0</span> <span class="ow">or</span> <span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">Q</span><span class="p">):</span>
            <span class="c1"># This outer plume segment is not viable</span>
            <span class="n">flag</span> <span class="o">=</span> <span class="kc">False</span>
            <span class="n">z0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">z_0</span><span class="p">,</span> <span class="n">z_lower</span><span class="p">])</span>
            <span class="n">y0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">len</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">len</span><span class="p">)])</span>
        
        <span class="k">else</span><span class="p">:</span>
            <span class="c1"># This outer plume segmet is viable...stop integrating widths</span>
            <span class="n">done</span> <span class="o">=</span> <span class="kc">True</span>
            
            <span class="c1"># Check where the diffuser is</span>
            <span class="k">if</span> <span class="n">z_lower</span> <span class="o">&gt;=</span> <span class="n">yi</span><span class="o">.</span><span class="n">z0</span><span class="p">:</span>
                <span class="c1"># This outer plume segment should not exist</span>
                <span class="n">flag</span> <span class="o">=</span> <span class="kc">False</span>
                <span class="n">z0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">z_0</span><span class="p">,</span> <span class="n">z_lower</span><span class="p">])</span>
                <span class="n">y0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">len</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">yo</span><span class="o">.</span><span class="n">len</span><span class="p">)])</span>
            
            <span class="k">else</span><span class="p">:</span>
                <span class="c1"># This is the next outer plume segment to integrate</span>
                <span class="n">flag</span> <span class="o">=</span> <span class="kc">True</span>
                <span class="n">z0</span> <span class="o">=</span> <span class="n">z_lower</span>
                <span class="n">y0</span> <span class="o">=</span> <span class="p">[]</span>
                <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Q</span><span class="p">)</span>
                <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Q</span> <span class="o">*</span> <span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">))</span>
                <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">s</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
                <span class="n">y0</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">p</span><span class="o">.</span><span class="n">rho_r</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="n">T</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
                <span class="n">y0</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">c</span> <span class="o">*</span> <span class="n">Q</span><span class="p">)</span>
    
    <span class="c1"># Return the results of the initial conditions search</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">z0</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">y0</span><span class="p">),</span> <span class="n">flag</span><span class="p">)</span></div>

<div class="viewcode-block" id="outer_fr"><a class="viewcode-back" href="../autodoc/spm/smp.outer_fr.html#smp.outer_fr">[docs]</a><span class="k">def</span> <span class="nf">outer_fr</span><span class="p">(</span><span class="n">u_0</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">bi</span><span class="p">,</span> <span class="n">rho_a</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">Fr_0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the outer plume initial width and velocity</span>
<span class="sd">    </span>
<span class="sd">    Computes the initial velocity of an outer plume segment using a Froude</span>
<span class="sd">    number condition analogous to Wueest et al. (1992) and calibrated and </span>
<span class="sd">    reported in Socolofsky et al. (2008).</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    u_0 : float</span>
<span class="sd">        Initial guess for the outer plume velocity (m/s).</span>
<span class="sd">    Q : float</span>
<span class="sd">        Flow rate in the outer plume (m^3/s).</span>
<span class="sd">    bi : float</span>
<span class="sd">        Radius of the inner plume (m).</span>
<span class="sd">    rho_a : float</span>
<span class="sd">        Density of the ambient water (kg/m^3).</span>
<span class="sd">    rho : float</span>
<span class="sd">        Density of the continuous outer plume fluid (kg/m^3).</span>
<span class="sd">    g : float</span>
<span class="sd">        Acceleration of gravity (m/s^2).</span>
<span class="sd">    Fr_0 : float</span>
<span class="sd">        Equilibrium plume Froude number for the outer plume (--).</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    u0 : float</span>
<span class="sd">        The initial velocity of the outer plume (m/s).</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># The Froude number condition is implicit; define the residual for use in</span>
    <span class="c1"># a root-finding algorithm</span>
    <span class="k">def</span> <span class="nf">residual</span><span class="p">(</span><span class="n">u</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Computes the residual of a modified Froude number condition for </span>
<span class="sd">        starting outer plume segments using the current guess for the initial</span>
<span class="sd">        velocity u.</span>
<span class="sd">        </span>
<span class="sd">        Parameters</span>
<span class="sd">        ----------</span>
<span class="sd">        u : float</span>
<span class="sd">            the current guess for the initial velocity (m/s).</span>
<span class="sd">        </span>
<span class="sd">        Notes</span>
<span class="sd">        -----</span>
<span class="sd">        All parameters of `outer_Fr` are global to this function since it is</span>
<span class="sd">        a subfunction of `outer_Fr`.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="c1"># Compute the outer plume width from Qo</span>
        <span class="n">b</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">Q</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">Q</span> <span class="o">*</span> <span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">bi</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        
        <span class="c1"># Calculate the deviation from the desired Froude number</span>
        <span class="k">return</span> <span class="n">u</span> <span class="o">-</span> <span class="n">Fr_0</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">((</span><span class="n">b</span> <span class="o">-</span> <span class="n">bi</span><span class="p">)</span> <span class="o">*</span> <span class="n">g</span> <span class="o">*</span> <span class="p">(</span><span class="n">rho_a</span> <span class="o">-</span> <span class="n">rho</span><span class="p">)</span> <span class="o">/</span> 
               <span class="n">rho</span><span class="p">))</span>
    
    <span class="k">return</span> <span class="n">fsolve</span><span class="p">(</span><span class="n">residual</span><span class="p">,</span> <span class="n">u_0</span><span class="p">)[</span><span class="mi">0</span><span class="p">]</span></div>

</pre></div>

          </div>
          
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
<h1 class="logo"><a href="../index.html">Texas A&M Oil spill / Outfall Calculator</a></h1>








<h3>Navigation</h3>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../user_manual.html">TAMOC User Manual</a></li>
<li class="toctree-l1"><a class="reference internal" href="../unit_tests.html">Unit Tests</a></li>
<li class="toctree-l1"><a class="reference internal" href="../glossary.html">Glossary</a></li>
</ul>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../bugs.html">Reporting Bugs</a></li>
<li class="toctree-l1"><a class="reference internal" href="../readme.html">Read Me File</a></li>
<li class="toctree-l1"><a class="reference internal" href="../release.html">Release Notes</a></li>
<li class="toctree-l1"><a class="reference internal" href="../license.html">License</a></li>
</ul>

<div class="relations">
<h3>Related Topics</h3>
<ul>
  <li><a href="../index.html">Documentation overview</a><ul>
  <li><a href="index.html">Module code</a><ul>
  </ul></li>
  </ul></li>
</ul>
</div>
<div id="searchbox" style="display: none" role="search">
  <h3 id="searchlabel">Quick search</h3>
    <div class="searchformwrapper">
    <form class="search" action="../search.html" method="get">
      <input type="text" name="q" aria-labelledby="searchlabel" />
      <input type="submit" value="Go" />
    </form>
    </div>
</div>
<script>$('#searchbox').show(0);</script>








        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="footer">
      &copy;2020, Scott A. Socolofsky.
      
      |
      Powered by <a href="http://sphinx-doc.org/">Sphinx 2.4.4</a>
      &amp; <a href="https://github.com/bitprophet/alabaster">Alabaster 0.7.12</a>
      
    </div>

    

    
  </body>
</html>